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Abstract. The main goal of this paper is to set up a numerical laboratory for the 
study of the slow evolution of the density and of the pressure tensor profiles of an 
j_j ■ otherwise coUisionless stellar system, as a result of the interactions with a minority 



component of heavier "particles". The effects that we would like to study are 
those attributed to slow coUisional relaxation and generically called "dynamical 
friction"; in real cases, or in numerical simulations, the processes involved are 
complex, so that the relaxation associated with the granularity in phase space is 
generally mixed with and masked by evolution resulting from lack of equilibrium 
or from a variety of instabilities and collective processes. We start by revisiting 
the problem of the sinking of a satellite inside an initially isotropic, non-rotating, 
spherical galaxy, which we follow by means of N-body simulations using about one 
million particles. We then consider a quasi-spherical problem, in which the satellite 
is fragmented into a set of many smaller masses with a spherically symmetric 
initial density distribution. In a wide set of experiments, designed in order to 
bring out effects genuinely associated with dynamical friction, we are able to 
demonstrate the slow evolution of the density profile and the development of 
a tangentially biased pressure in the underlying stellar system, while we briefiy 
address the issue of the circularization of orbits induced by dynamical friction on 
the population of fragments. The results of the simulations presented here and 
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others planned for future investigations allow us to study the basic mechanisms 
of slow relaxation in stellar systems and thus may be of general interest for a 
variety of problems, especially in the cosmological context. Here our experiments 
are conceived with the specific goal to clarify some mechanisms that may play a 
role in the evolution of an elliptical galaxy as a result of the interaction between 
the stars and a significant population of globular clusters or of the merging of a 
large number of small satellites. 

Key words, galaxies - dynamical friction - evolution of stellar systems 
1. Introduction 

If we consider a globular cluster of mass irigc ~ 10^ Mq orbiting inside a spherical galaxy, 
given the mass ratio 10^ with respect to the mass of a typical star, we find that it 
should suffer dynamical friction on a time-scale Tfr ^ lO^yr, through scattering of the 
stars. [The term dynamical friction refers to the slow relaxation process associated with 
discrete two-body encounters investigated in the pioneering studies of Chandrasekhar 
P943|l . In these classical analyses, estimate is given of the cumulative effect of discrete 
encounters on a test particle passing through a homogeneous system of field particles. 
This effect, associated with the granularity of the system in phase space, is ignored in 
pure Vlasov descriptions of stellar systems. The term dynamical friction is often extended 
to describe loosely the physical origin of orbital decays of satellites of finite mass also in 
real inhomogeneous systems and simulations, which fall obviously outside the classical 
description of Chandrasekhar.] If the cluster starts from a circular orbit located at r^, 
in a time ~ Tfr it will reach a lower orbit close to ry, while a fraction of the energy 
initially associated with such a cluster Ei — mgc^{ri) + {l/2)mgcri{d^ / dr)^ is released 
into the distribution function of the hosting galaxy. If we refer to the realistic possibility 
of the presence of several thousands of such globular clusters we may thus argue that 
on a fraction of the Hubble time the underlying stellar system may well absorb a non- 
negligible fraction of its binding energy and thus should slowly evolve (see also Lotz et 
al. lM?T|l . 

Under the suggestive scenario just outlined, the issue we would like to address in this 
paper is the following conceptually challenging problem. In a two-component system, 
where the more massive component (the "galaxy") is basically collisionless and the mi- 
nority component (the globular cluster system or a large number of mini-satellites caught 
in minor mergers) is slowly dragged in toward the galaxy center by processes of dynami- 
cal friction with the stars of the galaxy, how does the distribution function of the galaxy 
evolve during the process? A resolution of this challenging problem requires a clarifica- 
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tion of the basic mechanism of dynamical friction under inhomogeneous conditions and a 
discussion of the collective effects involved in the reaction of the underlying distribution 
function. Note that while the infall of a single cluster (or satellite) is a problem with a 
preferred direction (given by the angular momentum of the satellite), the capture of a 
large number of randomly distributed clusters (or mini-satellites) can preserve spherical 
symmetry, if present to begin with. 

It might be argued that the slow evolution in the quasi-isotropic many-cluster case 
described above could be modeled in terms of an adiabatic evolution of a distribution 
function within spherical symmetry. For the different problem of the adiabatic growth 
of a central massive black hole, such a procedure has indeed been developed (e.g., see 
CipoUina & Bertin [r994l and references therein). Here it is not clear how a semi- analytic 
procedure to calculate such slow evolution, taking advantage of the conservation of the 
relevant actions, can be formulated. In addition, given the relation between dynamical 
friction and resonant interactions (e.g., see Weinberg I1989|l . in spite of its slowness the 
process is likely to be inherently non-adiabatic. Therefore, it seems that the only viable 
way to an answer is that of approaching the problem by means of suitable N-body 
simulations. 

In the long run, results from this project could also shed light on another theoretically 
interesting scenario, which so far has found major applications only in the context of the 
internal dynamics of globular clusters. This is a picture, pioneered by Bonnor H1956|l 
in his analysis of gas spheres, by which stellar systems may be subject to a process of 
gravitational catastrophe (Lynden-Bell & Wood [n)68fl . when their central concentration 
exceeds a certain threshold value. The catastrophe consists in an instability process, 
where evolution leads to cores that become more and more concentrated. The source of 
this instability is thermodynamical, in the sense that it results from a counter-intuitive 
property of self-gravitating systems, which can be described in terms of negative specific 
heat. Recently, this general scenario has been demonstrated to hold for an astrophysically 
interesting family of partially relaxed stellar systems (Bertin & Trenti 2003). 

The study of relaxation in N-body simulations and by means of N-body simulations is 
a vast field of research (e.g., see Hohl llOfSl and Huang et al. 11993(1 . The study of the sink- 
ing of a single satellite inside a galaxy, as a result of dynamical friction, was undertaken 
systematically about two decades ago, initially through simplified simulations, within a 
restricted three-body problem and with the center of mass of the primary galaxy being 
kept fixed, and then by means of more realistic, self-consistent simulations (Bontckoe & 
van Albada lWzl Zaritsky & White UnHHI Bontekoe lH^ . The resuhs showed that the 
general scaling predicted by the classical formulae derived in the "homogeneous limit" (or 
"local limit") ( Chandrasekhar Fl 943|l was basically applicable, although the sinking time 
scale turned out to be significantly longer (by a factor of 2 — 3 for the regimes studied 
by the numerical experiments) if the self-consistent response of the galaxy was properly 
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incorporated. A convincing interpretation of the sinking process can be made, at least in 
the limit where the mass of the sinking satellite is small, by investigating the linear re- 
sponse in the galaxy distribution function through stellar dynamical techniques (Fridman 
& PolyachenkoEnHl Palmer & PapaloizouCnHSI Weinberg lTM)! [TM^ Bertin et al. HlMl 
hereafter BPRV94; Palmer ri994|l . Such response consists of an "in-phase" contribution, 
which has little relevance to the mechanism responsible for dynamical friction, and an 
"out-of-phase" contribution, associated with the "resonant" interaction, which produces 
a wake able to set a significant torque on the satellite. These techniques demonstrate 
that the collective wake in the galaxy, largely dominated by the dipolar (/ = 1) response, 
can trail significantly away from the orbiting object (when the satellite is outside the 
galaxy, the response is concentrated around a location opposite to the object with re- 
spect to the global center of mass), so that the resulting torque on the satellite that is 
responsible for dynamical friction can be significantly reduced with respect to other more 
naive calculations. Discrepancies between earlier simulations and the predictions of the 
stellar dynamic semi-analytic theory were ascribed to non-linearity effects (in the sense 
that the satellite-to-galaxy mass ratio Ms/Mq ~ 0.1 considered in the simulations was 
judged to be excessive; see Hernquist & Weinberg ll989|l . but the general picture appears 
to have been clarified. It is noted that the differences in behavior with respect to the 
naive expectations of the traditional formulae for dynamical friction become substantial 
when the size of the satellite is finite (Weinberg 1)1989(1 refers to the length ratio Rs/Rq', 
the multipoles associated with the perturbation with I above a certain threshold are 
bound to be unimportant), while the differences should become negligible for very com- 
pact satellites (in which case all multipoles contribute, with a divergence that is removed 
by the Coulomb logarithm). 

Several aspects of the general problem of dynamical friction have been revisited re- 
cently. The issue of the contrast between "local" and "non-local" effects has been re- 
discussed (Prugniel & Combes IBMI Maoz lBMl Cora, Muzzio & Vergne lTM?)! . with the 
confirmation that the description of the friction process in terms of the "homogeneous 
limit" fChandrasekhar ri943(l should be adequate in the limit of low-mass, relatively com- 
pact satellites. A very recent investigation (Cora, Vergne & Muzzio [200 l|l also suggests 
that the process depends very little on the regularity or the stochasticity of the stellar 
orbits in the galaxy, in contrast with some earlier conjectures fPfenniger I1986|l . Other 
interesting problems are the problem of circularization (of the satellite orbit), the rela- 
tion of the friction effect to the past history of the dynamical event under consideration, 
and the issue of a direct calculation of the relevant dissipative force from the response 
of the stellar system to the perturbing satellite (e.g., see Seguin & Dupraz [1994L 1199^ 
see also Colpi lDMI Colpi & Pallavicini lTMSI Colpi, Mayer & Governato lT^ Nelson & 
TremaineCnnnil- 
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Note that so far investigations of processes of dynamical friction have mostly focused 
on the issue of the fate of the sinking objects (and mostly for the case of one object; 
but see Tremaine, Ostriker & Soitzer ri975|l . while the problem of the reaction of the 
distribution function of the underlying stellar system (especially in the quasi-spherical 
limit of a system of many heavy objects outlined above) remains practically unexplored. 

Our paper is organized as follows. In Sect. [21 we summarize the theoretical framework 
that allows us to understand the mechanisms of dynamical friction within the context 
of a single sinking satellite. In Sect. Owe describe the model adopted. For the present 
exploratory investigation, the choice of the model (both for the galaxy and for the satel- 
lite or the fragments) is mostly dictated by convenience (for the study of dynamical 
mechanisms and for an easier comparison with previous analyses) rather than by real- 
istic requirements. In Sect. 21 we describe the diagnostic tools that will be used in the 
simulations and some preliminary tests. In Sect. [5l we proceed to illustrate the results of 
the simulations for the case of one sinking satellite and for the quasi-spherical case for 
which the sinking of many smaller mini-satellites is considered. In Sect.jHlwe briefly draw 
our conclusions and set out the plans for future investigations. 

2. Theoretical framework for the problem of a single sinking satellite 

Here we briefly recall the main issues involved in the problem of a single sinking satellite 
of "diameter" 2Rs, initially located on a circular orbit at distance sa Rq from the 
center of mass of the entire system. Here Rq denotes the "radius" of the galaxy. Let 
> be the angular velocity of the satellite along its circular orbit. The case of a very 
compact satellite is likely to be within the reach of the classical formulae for dynamical 
friction fChandrasekhar ri94H|l . In the case where the satellite is sufficiently extended, for 
the description of its density distribution and the associated gravitational potential, only 
the first harmonics in a multipole expansion are important, up to Imax ^ T^fs I {'^Rs) (see 
Weinberg CnHni). 

A "friction" torque acting on the satellite is expected to result from the action of 
suitable resonances with star orbits. To understand this, we may refer to the action of the 
satellite as a perturbation over a basic state that is non-rotating, spherically symmetric, 
and characterized by isotropic distribution function fo{E), where E — v"^ /2 + ^o{r) is the 
star energy, per unit mass, in the absence of the perturbation, and we may then follow 
the analysis developed earlier (BPRV94). 

Analytical results can be derived in the limit when the mass of the satellite is van- 
ishingly small (with respect to the mass of the galaxy, so that the resulting perturbation 
in the galaxy distribution function can be treated by a linear theory; see subsection 12. 2|1 . 
In practice, for real systems or for numerical simulations, some of the evolution effects 
are going to be associated with the finite mass of the satellite. Once we move to con- 
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sider the possibility of a quasi-equilibrium configuration for a composite system (galaxy 
plus satellite or a population of mini-satellites), we then have to worry about genuine 
instabilities that may occur, independent of the granularity present in the system. All of 
this will contribute to evolution although, in principle, it is not related to the original 
gedanken experiment at the basis of the classical definition of dynamical friction. 

Another subtle aspect of the problem is related to geometry. Broadly speaking, we may 
think that the difference between the case of straight orbits (considered in the classical 
description of dynamical friction) and that of orbits in the inhomogeneous potential of 
the galaxy is taken into account in the following discussion of the resonant response. Yet, 
the cumulative effects on a single star, because of multiple encounters with the satellite 
(for example, for stars on quasi-circular orbits), grow in time differently from the classical 
picture, even before dealing with the collective behavior of the entire set of stars that 
characterizes the resonant response. 

2.1. Single-star Resonances 

The frequency iv associated with the m-th component of a multipolar term of order I 
in the expansion of the potential of the satellite is readily related to the angular 
frequency of the satellite oj = —rnQg- Therefore, if we refer to Eq. (80) of BPRV94, we 
find that the appropriate resonance condition is —mfls—pflr{E, J) + sQ,o{E, J) = 0; here 
s, m, p are integers (— / < s,m < +1; — oo < p < oo) and flr{E, J) > 0, 0,e{E, J) > 
are the two frequencies (dependent on the star specific energy E and angular momentum 
J) that characterize the star orbit and depend on the properties of the basic state. As 
discussed in BPRV94, only the components with s = 1,1 — 2, .. — I contribute to the 
density perturbation. 

For an extended satellite, the p = component, which corresponds to an average 
along the star radial orbit excursion, is expected to dominate so that the resonance 
condition should reduce to 0,e{E, J) = rnQg/s. Note that in the course of time, because 
of dynamical friction, the angular frequency of the satellite is going to change. Clearly, 
the resonance conditions that we have just stated are applicable provided that 0^/0^ <C 
0,g{E, J), which we can rewrite as fig <C 0,^/1. 

2.2. "Shielding" and global resonances 

If wc consider the linearized Vlasov-Poisson equation for our problem, we sec that Eq. (54) 
of BPRV94 is modified into /i = C{^s + ^'i), where £ is a linear operator that basically 
describes one integration along the importurbed characteristics. The quantity /i can be 
called the linear response of the distribution function. This response has a "driven" part, 
associated with $„, and a "self-gravitating" part, related to the potential $i, that must be 
consistent with the response itself, following the Poisson equation V^$i = 47rG / fid^v. 
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In plasma physics, the latter contribution is often called the "shielding" associated with 
the collective behavior of the plasma. Such shielding is often ignored, although it should 
be emphasized that this step is improper, because the omitted term is a linear contribu- 
tion. A relatively simple analysis can then be made, by neglecting the self-gravity of the 
response, showing the connection between single-star resonances and dynamical friction. 
Note that this is the only ground where a comparison with Chandrasekhar's fl943l|l work 
can be made. 

The general case, which includes the possibility of resonance with discrete global 
modes of the system, is currently beyond the reach of analytical investigations (see also 
Vesperini & Weinberg 2000). 



3. Model 

As a model of the primary stellar system (the "galaxy"), following Bontekoe & van 
Albada (|1987|l we consider an isotropic polytropic model of index n = 3. The distribution 
function is thus of the form fair^v) = A{Ea - iJ)""^/^ for E < E^, with fG{r,v) = 
for E > Eq, where E = $g(»^) + t'^/S is the single-star energy, Eq = ^g{R) = -GM/R 
is the gravitational potential at the boundary r — R, and M is the total mass of the 
system. Polytropes can be described in terms of the Lane-Emden function 0{S,), which 
satisfies the differential equation: 

where n is the polytropic index, with the boundary conditions — I, d0/d^ = at ^ = 0. 
From the function 9{^) one can find the potential ^g{^) a-nd the associated density Pci''')- 

The secondary object (the "satellite") is described by a rigid Plummer density 
distribution (corresponding to a polytrope of index n — 5), characterized hy 9 ~ 
(1 -I- 0.5^^)~^/^. The potential and the mass distribution of the satellite are then defined 
by two scales, Rg and mass Ms, so that ^s{r) = -GMs{Rl -t-r^)"^/^, with cumulative 
mass Ms{r) — Msr^{R'^ + r2)^3/2^ rpj^^^ distribution extends to infinity, but 90% of the 
mass is contained within 3.7 Rg. 

As a result of the interaction with the satellite, the distribution function of the galaxy 
is a time-evolving / that will be associated with a mean field — V<i>, to be derived from the 
Poisson equation A(f>(r,t) — AttG J fd^v — 4'KGp{r,t). We solve the Poisson equation 
on a spherical polar mesh of 24 x 12 x 24 cells in the r, 0, and (p directions. The radial 
mesh is logarithmic. The angular mesh is fixed, with 12 equal divisions in cosfl and 24 in 
if. On this mesh the density, the potential, and the force field are expanded in associated 
Legendre functions P™(cos0) up to Pg. At each time-step the center of the spherical 
mesh used to calculate the potential and the force field in the galaxy is repositioned so 
as to coincide with the center of mass of the galaxy. The equations of the motion (for the 
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particles representative of the galaxy and, separately, for the satellite) are integrated in 
Cartesian coordinates with the so called leapfrog scheme. This is basically the N-body 
code described by van Albada & van Gorkom l|1977|l and then improved and used by van 
Albada H1H82(I and by Bontekoe & van Albada H1987|l . Therefore, we do not provide here 
a full description of the method employed, since it is available in the papers quoted. 



3.1. Sinking of a single satellite 

A first set of numerical experiments addresses the sinking of a single satellite, initially 
located on a circular orbit at 9 ~ 7r/2. 

The galaxy is considered as a collisionless stellar system, so that its representative 

particles (i — 1, , N) interact with each other through the mean gravitational potential 

$ and directly with the satellite: 

Note that <i> is updated separately by the Poisson-solver scheme, following the original 
method described by van Albada & van Gorkom H1977|l (this step effectively closes the 
system of equations), which guarantees a pure time- reversible evolution. 

In turn, the satellite is taken to interact with all the representative particles of the 
galaxy directly via two-body forces: 

d'^Vs _ Gm{r, - r^) 

In the equations above, r^, m are the position vectors and the mass of the representative 
particles of the galaxy, while Tg is the position vector of the satellite. All position vectors 
are measured with respect to an inertial frame of reference with its origin at the center 
of mass of the combined (galaxy -I- satellite) system. 



3.2. Sinking of many fragments 

Another set of experiments studies an initial configuration where the satellite is frag- 
mented into several {Nj) pieces ("mini-satellites" or "fragments"). For simplicity, we 
take the fragments to be characterized by the same Plummer density distribution with 
mass Mf (such that Ms — NjMf) and the same lengthscale Rf = Rs (a smaller length- 
scale Rf would be more realistic for the astronomical system we have in mind, but would 
be less easy to simulate). 

The mutual interactions between mini-satellites can be kept "on" by considering 

dt^ ^ ^ ^' ^ ^ [i?^ + (r.-r/,)2]3/2' 
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for the representative particles of the galaxy (labeled by i = 1, iV), and 



(5) 



for the mini-satellites (labeled hy j = l,...,Nf), where 



i=N 



a 



(1) ^ 



(6) 




is the acceleration produced by the galaxy and denotes the gravitational acceleration 
exerted on the j-th niini-satcUitc by the other fragments. The gravitational interaction 
between two extended bodies is complicated to describe. For the purposes of the present 
paper, where we are not interested in the phenomena associated with the details of such 
interaction, we simplify the simulations by referring to the following prescription; 



All position vectors are measured with respect to an inertial frame with its origin at the 
center of mass of the entire system. 

The mutual interaction between mini-satellites can be ignored and turned "off' simply 
by taking <x?^ = 0. This procedure may be the more realistic procedure if we wish to 
simulate a system of globular clusters inside a galaxy. In fact, for such system, the few- 
body-problem effects associated with such interaction are expected to play a secondary 
role. 

As to the initial distribution of the fragments we have considered three cases, to be 
described below. The first is the most natural generalization of the study of the sinking 
of a single satellite. A posteriori, we have found that in such a case the effects related to 
the lack of equilibrium in the initial conditions are too strong and confuse the picture of 
evolution. (Similar problems, related to the lack of equilibrium, were present also in the 
case of a single satellite, but those are not emphasized in this paper because that case 
has been well studied in the past and we prefer to focus here on the new quasi-spherical 
problem in the presence of many fragments). Therefore we have devised quieter starts 
so as to better disentangle the effects associated with the granularity of the system from 
those associated with the lack of equilibrium. 

Given the relatively small number of fragments involved, we anticipate that differ- 
ent realizations of the physical cases described below may be associated with different 
evolutions. 




(7) 



3.2.1. I: Fragments quasi-isotropically distributed on circular orbits on a thin shell 

We first consider the case where all the fragments are distributed initially with random 
positions on a thin shell, located at radius rsft,eH(0)- The fragments are given initial veloc- 
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ities appropriate for circular orbits, with random orientations, based on the gravitational 
field generated by the galaxy alone. 

The purpose of simulations of this type is to derive information relevant for a more 
realistic situation where the fragments have a full three-dimensional space distribution, so 
that the interaction between the galaxy and the minority population is differential with 
radius. Such ideally desired simulations are not feasible with current codes. To move one 
step further in the desired direction we have then considered the case of two shells of 
fragments initially located at different radii. 

Unfortunately, due to the finite mass of the set of fragments, these simulations are 
characterized by violent epicyclic oscillations that shake the system on the fast dynamical 
time scale and tend to persist long after the initially expected transient. 

3.2.2. II: Fragments quasi-isotropically distributed on circular orbits on a thick shell 
with improved initial conditions 

In order to set up a smoother quasi-equilibrium initial condition, we have decided to 
proceed as follows. We refer to a spherically symmetric density distribution for a shell of 
the form psheii = Poexp[-(r - r shen{^)Y / a?] for |r - rsheu{^)\ < 2a (and vanishing for 
— fsheiiiO)\ > 2a). The normalization constant po is taken in such a way that the total 
mass associated with the shell is Ms- In view of the choice of model for the fragments, 
we have decided to take 2a = Rs- 

Such density distribution generates a potential ^sheiii'''), which can be calculated 
numerically. We now consider a modified distribution function for the galaxy fQ{r,v) 
taken to be of the polytropic form as in Sect.OJ but with E = $g(?') + ^sheiiir) + jl. 
We then have to integrate a Lane-Emden equation similar to Eq. but modified by the 
presence of the external potential ^sheiiif)- This yields a potential ^gI?") different from 
that considered in Sect.|31 to be inserted into the definition of /g(^, v). 

Finally, we consider the clumpy realization of the shell density distribution by reducing 
it to a set of N f fragments of mass Mf. The j— th fragment is initially distributed in r at 
a location r fj (with random angular coordinates) so that the set of fragments reproduces 
the assumed density distribution of the shell; its initial velocity is that appropriate for 
the circular velocity of a test particle in the combined potential ^cii^) + ^sheii{r)- 

This procedure ensures a quieter start while leaving the basic picture of the single 
shell case unchanged. It can be generalized to the case of two shells of fragments initially 
located at different radii. 
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3.2.3. Ill: Fragments constructed by clumping part of the galaxy distribution 
function 

The cases considered above are aimed at describing the interaction between the galaxy 
and a minority component with phase-space properties distinct from those of the galaxy. 
As we have seen while addressing the need for a quieter start, such distinction is likely 
to be accompanied by effects that are due to the lack of equilibrium in the initial condi- 
tions, difficult to separate from the secular effects that are traditionally called dynamical 
friction, fn this assessment, we can go even one step further and we may argue that 
some of these initial conditions (with the different phase-space properties considered) 
may actually be characterized by dynamical instabilities. If this were the case, it would 
be very hard or even impossible to distill the true effects of dynamical friction out of 
the simulations. In order to address this point, we will run parallel simulations of case 
II with a smooth, Vlasov realization of the shell density, which should allow us to single 
out the effects of possible Vlasov instabilities of such system (see Subsect. 13.2. 4|) . 

In view of this discussion, we have decided to address a third set of simulations that 
is less relevant for the astrophysical case that has motivated this paper but is definitely 
more interesting for the study of the problem of dynamical friction. 

Here we consider the initial equilibrium state to be basically that described by the 
galaxy model alone, without satellites or shells or fragments, as defined in the first para- 
graph of Sect. 13 The mini-satellites here are then introduced by clumping together part 
of the distribution function fa in such a way that it corresponds to a diffuse "shell" in 
space. Formally, we separate out a piece of the distribution function f shell (for simplicity, 
we keep the notation "shell" even if we are considering a situation different from that 
of case I and case II). Then the system made of /g — f shell is simulated by the Vlasov 
code, while the contribution fsheii is reduced to clumps that are brought back to the 
form of the fragments and treated by means of direct interaction with the galaxy repre- 
sentative particles (as described in subsection 13. 2|l . Note that in this case the "clumps" 
or "fragments" are distributed in phase space exactly as the particles of the galaxy (thus 
minimizing the lack of equilibrium in the initial conditions and the risk of introducing 
undesired sources of instability); in particular, the clumps are initially characterized by 
an isotropic distribution in velocity space. 

We should emphasize that this third class of simulations, while addressing the settling 
of heavy masses dragged inwards by dynamical friction, is conceptually distinct from 
the study of the mass stratification process in self-gravitating systems, which may be 
performed by means of direct or Fokkcr-Planck codes (e.g., see SDitzer [l987|l . 
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3.2.4. Purely collisionless reference models 

For the three cases just outlined, the simulations carried out following the description 
given in subsection 13 . 21 will be accompanied by reference simulations of similar but purely 
collisionless models. These can be seen as cases for which the number of fragments for- 
mally diverges, while the total mass contained in the fragments remains constant. In the 
simulations the case of very large number of fragments cannot be treated in terms of di- 
rect interaction; instead, for these cases, the fragments are treated as simulation particles 
of the collisionless component. It is clear that for case III, in such reference collisionless 
models, the identification of f shell is meaningless. In practice, to monitor the properties 
of our integration scheme, we have labeled the representative particles associated with 
f shell by a different "color" . 

The purpose of these parallel reference runs is to check, as much as possible, that we 
are not confusing evolution effects associated with the initial conditions with the effects 
that are genuinely associated with the granularity of the system. 

3.2.5. Initial set up 

In order to derive the initial distribution of the particles in the simulation we start from 
the mass distribution M{r) which is known on a discrete set of points labelled by the 
index j, with Mj = M{r^). 

First we determine the distance of a particle from the galaxy center. If all the particles 
of the galaxy have the same mass m ~ M/N, then for the i*^ particle we define xx — 
m ■ {i — 1/2), and then find the cell j such that Mj > xx > Mj_i. The distance from 
the center of the i^^ particle is then found by an interpolation between r^^^,r^, and 
r^^^ . Next we choose two uniformly distributed random numbers: q,^ € [0, 1] for the cp 
coordinate of a particle, such that tp = 2nqi^, and qg G [—1,1] for the 9 coordinate, such 
that cos 6 — qe. The Cartesian coordinates are then found from the spherical coordinates 
by standard formulae. The choice of the initial velocities is performed in a similar way, 
using a rejection method (e.g., see Press et. al. ()1992|l ') to obtain the desired distribution 
function. For the potential <i>G(?') we proceed as for the mass distribution; the value of 
the potential at the radius of the particle location is found by interpolation between the 
values at r^~^,r^ , and r^^^ . 

In Case III, of a "subtracted" shell, the particles of the galaxy are divided in two 
"species": light particles, with mass m — M/N — {M — Ms)/Ni, and heavy particles 
with mass M f = Mg/Nf respectively, where Ms is the mass of the shell, Nf the number 
of heavy particles, and Ni the number of light particles. Initially, the heavy particles are 
all located inside a thin shell. Inside this shell there are no light particles. The i*^ light 
particle, for i e [l,il = iV • Al (r shell) /M], where rsheii is the position of the "shell", is 
located as in the non-subtracted case at the distance obtained by interpolation between 
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the r^~^,r^, and r^+^, where j is such that the condition Mj > m ■ [i — 1/2) > -Aij-i 
is satisfied. Then, for the /c"* heavy particle, with k E [l,Nf], we adopt xx — m ■ (il — 
1/2) + Mf {k — 1/2) and then find the cell j where Alj > xx > Mj-i. The distance fi-om 
the center of the k^^ heavy particle is then found by the interpolation rule mentioned 
above. After that the remaining light particles are distributed according the previous 
rule. Finally, in the simulations, the light particles interact with one another through 
their mean field, while the heavy particles interact directly with the light particles and 
among themselves, as described at the beginning of Sect. 13.21 

3.2.6. Other possible initial conditions 

Another possible initial condition is to turn on slowly the gravity associated with the 
massive objects or massive shells. In a collisionless simulation, it would lead to the "em- 
pirical construction" of equilibrium systems in which dynamical friction is absent and 
would thus be, for our purposes, essentially equivalent to the procedure described in 
Case II. Such construction of equilibrium models can also be approached analytically 
(see Cipollina & Bertin [T994I and references therein). 

The case where the satellite or the fragments are treated as discrete objects would 
provide information on the friction force, different from the information provided in this 
paper, but such approach would have no special merits, in terms of physical justification 
and physical interpretation, with respect to the one adopted here. 

3.3. The choice of the code 

The choice of code that we have made (basically, the code used by Bontekoe & van Albada 
appears to be appropriate for the study of the effects of dynamical friction, which 
are dominated by low-Z wakes (see Sect. 2). Still, one may ask the reason why we have 
resorted to that code instead of codes that are commonly utilized for a variety of studies 
in galactic dynamics or in the cosmological context. In particular, one popular class of 
codes, the tree-codes (e.g., see Barnes & Hut ()198f)|l \ has found modern realizations 
(e.g., GADGET; see Springel, Yoshida, & White I7()()l|l that are particularly appealing 
since they are flexible, widely used, and have gone through many tests. Tree-codes treat 
near particles in terms of direct interactions via a softened potential, because the use of 
the true Newtonian potential would introduce a generally undesired coUisionality among 
the simulation particles. Tree codes are often used to study the evolution of collisionless 
systems; for this purpose, the softening length has to be chosen in an optimal way, 
depending on the physical characteristics of the system that is being investigated. 

The use of a tree code in our project would then pose severe interpretation problems, 
because the non-physical effects associated with the introduction of the softening param- 
eter would show up precisely in relation to the issues that we are trying to investigate. 
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4. Diagnostics and test cases 

The results of our simulations can be studied in terms of the energy and angular mo- 
mentum budget as a function of time and by means of figures that illustrate the spatial 
structure of the density and of the pressure distribution within the galaxy in the course of 
time. Another set of figures will provide plots of the azimuthally averaged density profile 
for the galaxy and of the location of the satellite as a function of time; for the case 
of many mini-satellites, we will show a plot of the radius of the sphere that contains 
different fractions [x percent) of the total mass in the form of mini-satellites. We first 
proceed by specifying units and definitions for the relevant quantities. 

4.1. Units 

Unless specified otherwise, the units adopted for mass, length, and time are the mass 
of the galaxy (M), the radius of the galaxy (i?), and the natural crossing time {tcr — 
GM^/'^{2Kgai)'^/^ evaluated at the be ginning of the simulation; here Kgai is the total 
kinetic energy associated with the stars in the galaxy). For the galaxy model considered 
in this paper, the revolution period relative to a circular orbit at the periphery (at R) is 
« 11 tcr- To return to physical units, we may refer to the case where M = 2 • 10^^ Mq, 
R = 10 kpc, so that tcr — 0.18 • 10^ yrs and the revolution period at the periphery is 
« 2 • 10^ yrs. 

4.2. Conservation laws 

In the course of the simulation the total energy and the total angular momentum should 
remain constant at their initial values Etot and Jtot- 

In the case of a single satellite we write the energy conservation as 

Etot — Egal + Ksat + Eint = (8) 

1=1 

where Egai — Kgai + Wgai is the total (kinetic plus gravitational) energy of the galaxy 
in the absence of the satellite. From the point of view of the satellite, it is natural to 
introduce the energy Esat = Ksat + Eint, so that we may split the total energy into 
Etot — Egal + Esat- The angular momentum conservation is given simply by Jtot = 

J gal + J sat = Ei=f ™ i ^ ""i + ^s Vs X Vs . 

In the case of several mini-satellites, following the model outlined in Sect. |31 the 
energy takes the form 

i=JV 

£^404 = Egal + {Kfra + Wfra) + E.^t = ^ ^ ^'"^ ~^ ^^'''^^ ^ 
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while the angular momentum can be written as J tot = Jgai + J/ra = Si=i^ mri x Vi + 
X^j=i'^ '"'fj ^ Jj ■ Note that the energy associated with the system of mini-satellites 
now includes a term Wfra that describes the gravitational energy internal to the system 
of fragments. By analogy with the case of the single satellite, we may wish to refer to the 
energy Ejra = [Kjra + Wfra) + Eint as the total energy associated with the fragments, 

so that Etot — Egal + Efra ■ 



4.3. Density and pressure distributions 

In order to describe the galaxy density perturbation, we refer to the density response 
defined as pi{t, r) = p{t, r — ro{t)) — po{r); here po{r) is the adopted density distribution 
of the galaxy in the absence of satellites, ro{t) is the position of the center of mass of 
the galaxy in an inertial frame of reference. To diagnose the structure of the response it 
will be useful to consider the characteristics of the first multipoles. 

In the simulations the density on the spherical mesh is expanded in associated 
Legendre functions P™(cos^) up to Z = 6 

1—6 r- rn—l 

p{r,0,^) = ^ Ai,o{r)Y,^{cosO) + ^ V2A^,^(r) Re(y,-(cos ^, ^)) (10) 



1=0 



m=l 



m=l 



+ J2 V2Bi^m{r)lm{Yr{cose,^)) 



m—l 



where the coefficients of the expansion into real spherical harmonics are given by suitable 
integrations over the relevant solid angle 



(11) 
(12) 

(13) 



Ai,m{r)=V2 J p{r,e,ip)Re{Y{^ {cos 0,ip))dn , 

Aiflir) = J p{r,e,<p)Yi° {cos e)dn , 

Bi,m{r) = p{r,e,ip)Im{Yi"'{cos9,ip))dn . 

In particular, the dipole {I = 1) component of the density response pi is given by 

pl=i{r,0,ip) = Alfi^/3/Ancos6 - Ai_i-\/3/47rsin^cos<^ - y/S/An sin 6 sin (p . (14) 

The galaxy pressure tensor is defined in terms of the star distribution hmction in the 
standard way. In the case of a single satellite sinking into a galaxy along the equatorial 
plane we have {vr) — and {v0) ~ 0. 
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4.4. Tests 

We have performed many preliminary experiments aimed at testing the numerical code. 

The numerical implementation of the polytropic equilibrium configuration of Sect.|31 
with N — 5 ■ 10'^ ~ 5 • lO'^ particles has been tested by generating an unperturbed galaxy 
that has remained quiescent for about 50 crossing times. The radii of the spheres con- 
taining 0.5%, 1%, 1.5%, ....,5% of the total mass have been checked to remain constant 
in time to within 1%. 

Then, we have tested that a satellite with Rg = 0.1 R and very small mass, Alg = 
10~* AI, behaves as a light "test particle" without disturbing the main galaxy. This 
satellite, placed initially at the periphery of the galaxy = i? on an orbit with Vs,p — 
0.536 R/tcr and Vgr — 0.954 • 10^"^ R/tcr, remains very close to its initial circular path, 
never exceeding an eccentricity of e = 0.05. At the end of the run, at time t ^ 50 t^r the 
orbital energy of this light satellite is conserved to within 1%. 

4.4.1. Energy and momentum conservation 

In addition to the specific tests mentioned above, we have checked that globally energy 
and angular momentum are conserved, at t ~ 50 tcr, to better than 1%. 

4.4.2. Comparison with earlier work by Bontekoe and van Albada (1987) 

The main results of the simulations of the sinking of a single satellite by Bontekoe & 
van Albada (1987) have been checked in detail, with the use of simulations with N up to 
5 • lO'^ (see Sect. 15 .ll below). This adds much confidence, in relation both to the overall 
structure of the code used in this paper and to the modeling of dynamical friction used 
in these investigations (since results have been checked to be largely independent of N 
over a range significantly larger than that available earlier). 

5. Results of the simulations 

We have performed a number of runs, mostly based on iV = 5 • 10^. Runs labeled by A 
refer to experiments with a single sinking satellite, initially located on a circular orbit at 
rg = R; various combinations of satellite mass Mg/M and satellite size Rg/R have been 
considered. For these runs (some are listed in TableCJ, we have recorded the time tfaii /t^r 
at which the satellite reaches the central regions; for some experiments (in particular, 
for those labeled As and Ag) the satellite is unable to reach the center and we have thus 
studied the radius of minimum radial approach r,„i„. 

Runs labeled by B are listed in Table [21 and refer to the case of the slow evolution of 
the system in the presence of many fragments, initially located in a single shell, following 
the procedure of case I described in Sect. 13.2.11 Runs of type BS are based on the 
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smoother initialization of case II outlined in Sect. 13.2.21 Runs of type BT refer to the 
procedure of case III of Sect. 13.2.31 Here the fourth column records the initial radius of 
the shell and the fifth column lists the time at the end of the experiment; all runs are 
made with the mutual interaction, among fragments, on, but we have also performed a 
few runs for which the interaction is turned off. Runs of type C involve two shells of 
fragments. 



Run 


Ms/M 


Rs/R 


tfall/tcr 


Notes 




0.1 


0.1 


48.125 




Ae 


0.05 


0.1 


110.0 




As 


0.05 


0.05 


73.15 


Tmin = 0.05 R 


As 


0.05 


0.025 


60.5 


Trmn = 0.022 R 



Table 1. Single satellite, rs(0) = R. 



Run Nf Mf/M rsheii{0)/R tfin/tcr Notes 



Bi 


100 


0.001 


1 


220 


case I 


B2 


100 


0.001 


0.2 


110 


case I 


B3 


100 


0.001 


0.4 


110 


case I 


Bi 


100 


5 • lO""* 


1 


220 


case I 


Bs 


100 


5 • 10-* 


0.2 


110 


case I 


Be 


100 


5 • 10-* 


0.4 


110 


case I 


BSi 


20 


0.005 


0.3 


220 


case II 


BS2 


100 


0.001 


0.3 


220 


case II 


BSa 


20 


0.005 


1. 


110 


case II 


BSi 


100 


0.001 


1. 


110 


case II 


BTi 


20 


0.005 


0.3 


220 


case III 


BT2 


100 


0.001 


0.3 


220 


case III 


BT3 


20 


0.005 


0.7 


110 


case III 


BTi 


100 


0.001 


0.7 


110 


case III 


Ci 


200 


5 • 10-* 


1, 0.2 


220 


case I 


C2 


200 


5 ■ 10-* 


1, 0.4 


220 


case I 



Table 2. Simulations with shells of fragments. 
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5.1. The case of one sinking satellite 

We have studied the energy transfer between the galaxy and the sinking satelhte and the 
overall energy balance, using the total energy conservation as a diagnostic of the code 
accuracy. In particular we have considered the quantity AEgat — Ksat + Eint — Ksat (t — 
0) — Eint{t — 0). For example, for run Ai (characterized by M^/M = 0.1, Rg/R = 0.1, 
tfaii/tcr — 48.125) the energy lost by the satellite is initially very small (on the order of 
0.02, in the units following the conventions described in Sect. 14. ll at time t w 38.5 tcr), 
and then rapidly rises (to 0.075 at time t « 49.5 <cr); during the run the total energy 
is conserved to better than 0.002. Similarly, we have considered the angular momentum 
balance. Here the relatively large amount of angular momentum already associated with 
the galaxy right after the beginning of the simulation {Jgai ~ 0.0054 at i « 5.5 tcr) is 
not yet due to the scattering of star orbits by the sinking satellite, but is mostly related 
to the overall orbital motion of the galaxy in the inertial frame of reference of the center 
of mass; eventually, the satellite is dragged in, down to the center, and loses completely 
its angular momentum, which is acquired by the galaxy. 

During the process, following the steps indicated by Bontekoe & van Albada (1987), 
we can "measure" the Coulomb logarithm. The measurement of In A consists of the 
following steps: (i) Determine the velocity Vs of the satellite with respect to the local 
streaming motion of the stars (the streaming motion is determined as a mean veloc- 
ity of the stars in the "vicinity" of the satellite); (ii) Measure the energy Egat of the 
satellite at many instants, then determine dEgat/dt; (iii) Determine the density pc of 
the galaxy at the position of the satellite by linear interpolation of densities between 
the 8 nearest mesh corners; (iv) Calculate F{vg), the fraction of particles with ve- 
locity smaller than Vg. Finally the Coulomb logarithm is estimated from the relation 
47rG2 In A = -Vg{dEgat/dt)/{M^pGF{vg)). 

The main results that we have obtained by examining the experiments for the case 
of a single sinking satellite are the following: 

— Qualitatively, the general conclusions derived by Bontekoe & van Albada (1987) have 
been checked to hold. However, we have found that, while the classical formula of 
dynamical friction can be used to fit the orbital decay of the satellite over a significant 
time interval, the determined value of the Coulomb logarithm changes significantly in 
the course of the simulation (see Fig.^; during the initial transient and the final part 
of the orbital evolution, when the satellite reaches the central regions, the change in 
the nominal value of the Coulomb logarithm can be dramatic, as already noticed by 
Bontekoe & van Albada (1987). 

— The process of dynamical friction has been checked to operate via the generation of 
wakes, of which the dipolar I = 1 component appears to play a dominant role (see 
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Figs. 121 and 13). This generally agrees with the scenario outlined by Weinberg (1989); 
see also Sect. |51 

— As a result of the scattering of stars due to the sinking satellite, the galaxy acquires 
not only some angular momentum (the induced systematic motions are approximately 
those of solid body rotation, see Fig.^, but also a significantly anisotropic pressure 
tensor. Figure relative to run Ai, suggests that the final configuration reached 
is reasonably well characterized by Prr ~ P09, while the toroidal pressure p^^ is 
significantly larger even in the central regions. As a result, except for the outer parts, 
the evolved distribution function of the galaxy may turn out to be accessible to models 

with / = /(£;, J,). 

— In some runs, for which the diameter of the satellite is reduced, the satellite turns out 
to fall inwards faster, but to be unable to reach the center of the galaxy (see Fig.O. We 
have checked that this effect does not depend on the number N by running a simulation 
of type Aq with N up to 2 • 10^ (in such a way that the number of simulation particles 
in the region inside the radius of minimum approach increases from about 600 to 
about 2500; correspondingly, the observed radius of minimum approach changes by 
less than 1 %). We should recall that simulations with a code of the type we are 
using, with a finite expansion in spherical harmonics, cannot deal realistically with 
spatially small satellites (see comment at the end of Sect. 12. l|) . Still, a "saturation" of 
dynamical friction of the kind noted in our simulations is likely to be real, because the 
modeling problem associated with the spherical harmonic expansion is less severe for 
the late phases of the simulation, when the satellite is close to the galaxy center. The 
ultimate fate of a single sinking object might be related to other interesting issues, 
such as those met in the discussion of the formation of massive black hole binaries 
(e.g., see Quinlan & Hernauist 119971 Milosavljevich & Merritt [200 l|l . which go well 
beyond the scope of the present paper (that addresses the slow evolution of a galaxy 
on the large scale as a result of the capture of a large number of small objects); 
what we have noted here is just meant to say that on this specific feature we confirm 
the results of previous investigations (in particular, those of Bontekoe & van Albada 
1987) carried out with a much smaller number of particles. 

— Finally we note that the overall density distribution of the galaxy, at the end of the 
simulation, turns out to be less concentrated with respect to the initial distribution 
(see Fig.[7||. In passing, we note that a figure such as Fig. [T] naively suggests that the 
total mass is not conserved, but we have checked that in this respect this is only a 
false impression left by our way of representing the density profile; in particular, one 
should recall that the weight of the volume element rapidly increases with radius. 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 
M(r)/M M(r)/M 



Fig. 1. Numerical measurement of the Coulomb logarithm (left) and of the satellite energy 
loss by dynamical friction (right) for run Ai (characterized by Ms/M = 0.1, Rs/R = 0.1, 
tfaii/tcr = 48.125). The quantities are plotted against a Lagrangian radial coordinate in the 
same format as in the paper by Bontekoe & van Albada (1987). 

5.2. The case of many fragments 
5.2.1. Case I 

For various runs of type B and C we liave cliecked tliat the total energy conservation 
generally holds, over a period of « 100 tcr, to better than one part over 500. Yet, it is 
difficult to estimate the impact of dynamical friction on the galaxy evolution, because 
secular effects are masked by effects related to the initial lack of equilibrium. Prom this 
class of numerical experiments we have obtained the following results: 

— Wc have observed a process of shell diffusion, which is a collective effect qualitatively 
different from what might have been imagined by simply superposing the effects 
of orbital decay of the individual fragments. The phenomenon persists, presumably 
via interaction with the wakes induced by the fragments in the galaxy, even when 
the direct fragment-fragment interaction is turned off. Sharp coherent oscillations 
are noted at early times; they are recognized as epicyclic oscillations generated by 
our choice of initial velocities for the fragments. We have also performed tests on 
the process of shell diffusion and oscillations on the basis of 400 fragments initially 
located close to r = 0.6 R. 

— We have noted systematic changes of the galaxy density distribution. The general 
trend emerging from the simulations is that significant peaking of the density distri- 
bution may occur, provided the initial radius of one shell of fragments is well inside 
the galaxy. However, much of the change is apparently associated with the initial lack 
of equilibrium. 

— The process of density profile modification is not accompanied by the development of 
significant pressure anisotropy. 
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Fig. 2. Six frames describing the equatorial density response p\ = p — po oi the galaxy to 
the infalling satellite for run A^. The density wakes are portrayed at time t/tcr = (a) 38.5, (b) 
39.875, (c) 41.25, (d) 42.625, (e) 44, and (f) 53.625. The white circle denotes the location of the 
satellite. 



In conclusion, simulations of Case I turn out to be of little direct use for the main 
objectives of our project. However, they are extremely useful and instructive, because 
they demonstrate that, in order to properly single out slow relaxation processes, one has 
to be very careful about the choice of initial conditions and, in particular, because they 
provide concrete examples of some misleading effects that can occur as a result of the 
use of improper models. 
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Fig. 3. For each frame of Fig.|21we illustrate in detail the monopole {I — 0) and the dipole {I = 1) 
contributions to the density wake in the galaxy (the various curves represent the coefficients Aio, 
^11, ^00, and Bn defined in the text; see Sect. l4.3(l . 



5.2.2. Case II 

The use of the smoother initiahzation described as case II in Sect. 13.2"^ leads to runs 
where the secular effects associated with the granularity of the fragment population 
are more convincingly identified. The orbital decay of the fragments is illustrated in 
Fig. |S1 where we also show for comparison a parallel run where the shell is treated as 
collisionless. The interaction which generates dynamical friction of the galaxy on the 
fragments is responsible for the slow, but systematic, development of a tangential bias 
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Fig. 4. For run Ai we plot the mean azimuthal velocity of the galaxy measured on the equatorial 
plane at three different times; the rotation induced by the angular momentum exchange with 
the infalling satellite is approximately that of a solid body. As usual, time is measured in units 

of tar- 
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Fig. 5. For run A\ we plot two othogonal equatorial cuts of the anisotropy distribution generated 
in the galaxy by the infalling satellite at time t = 53.625 tcr, after the satellite has reached the 
center. 



in the pressure tensor, as illustrated in Fig. |3 and for an overall change in the density 
profile of the galaxy (see Fig. I1(J|) . 

The check against the coUisionless run is worth a special digression. In order to di- 
agnose the collisionality present in our system better, we have performed the following 
test. We have studied the correlation between single particle energy (for the particles 
simulating the galaxy) at two different times. If the simulation is truly coUisionless, 
the single particle energy Ep = mv^ jl + m<i>(r), where ^(r) is the mean field po- 
tential, should be strictly conserved. Indeed, by means of a plot similar to the lower 
right frame of Fig. |H1 we have checked that for the run illustrated by the lower left 
frame of Fig. |S1 more than 90% of the particles have an energy correlation such that 
\\Ep{t = 118 tcr) - Ep{t = 14 t^r)\lEp{t = 14 tcr)| < 0.06. In turn, for run BS2, oth- 
erwise similar, but characterized by the presence of 100 fragments, a similar level of 
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Fig. 6. The curves illustrate the orbital decay process for the satellite for three different runs 
{Ae, As, and Ag, labeled as 1, 2, and 3 respectively) characterized by different spatial size of 
the infalhng satellite {Rs = 0.1, 0.05, 0.025 R; Ms/M = 0.05). 




-0.4 -0.2 0.0 0.2 0.4 X 20 40 60 80 t/t 



Fig. 7. Left: The density distribution change induced in the galaxy by the infall of the satellite 
(run ^i). The thick upper curve (1) represents the initial distribution, while the thin lower curve 
(2) is associated with the final state. Right: Orbital decay for run Ai (cf. Fig. |5J as in Fig. |S| 

correlation for single particle energy (at 6 %) is preserved within the same time interval 
only by 50 % of the particles; the lower right frame of Fig.|Hlthus illustrates the "damage" 
associated with dynamical friction. 

5.2.3. Case III 

The even smoother procedure of case III leads to slower orbital decays, as illustrated 
in the lower frames of Fig. 1111 For this type of simulations, the initial distribution of 
fragment orbits replicates the isotropic star distribution of the underlying galaxy. Because 
of dynamical friction, the fragment orbits are subject to a slow process of circularization, 
as illustrated in the upper frames of Fig. 1111 

For each of the fragments that are involved in the simulation described as Case II, we 
can repeat the analysis that has led to the measurement of the Coulomb logarithm from 
the study of the decay of the orbit of a single satellite. For example, for run BS2 we have 
100 independent such determinations. In Fig. 1121 we record, for run BS2 a measurement 
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Fig. 8. The fall of a shell of fragments towards the galaxy center, for the case of Nf = 20 
(left frame on the top, run BSi) and for Nf — 100 (right frame on the top, run BS2)- For 
comparison, the left frame at the bottom describes the time evolution of a coUisionless shell (a 
run with Nf = 5 ■ 10*). Displayed lines represent 5, 15, 25, 35, 45, 55, 65, 75, 85, and 95 % of 
the shell mass. The right frame at the bottom displays, for run BS2, the correlation between 
single particle energy at time t — 14 tcr and energy at time t — 118 tcr, showing the effects of 
scattering associated with dynamical friction; units for energy are those indicated in Sect. 14. ll 
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Fig. 9. The development of pressure anisotropy in the galaxy as a result of the interaction with 
a shell of fragments dragged in towards the galaxy center, for the case of Nf — 20 (left frame, 
run BSi) and for Nf — 100 (right frame, run BS'2). The dashed curves represent the evolving 
value of Kt/2, where Kt is the total kinetic energy associated with the star motions in the 
tangential directions; the solid curve represents the evolving value of Kr, the total kinetic energy 
associated with the star motions in the radial direction. 
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Fig. 10. The change in density profile for the galaxy as a result of the interaction with a shell 
of fragments dragged in towards the galaxy center, for the case of Nf — 20 (run BSi). The 
thicker curve (1) represents the initial density profile, the thinner (2) gives the final profile. 
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Fig. 11. The effects of dynamical friction on a shell of fragments created out of the initial 
distribution function (case III of Sect.j^l- The lower frames show the general decay of orbits of 
the fragments for the case Nf = 20 (run BT3) and Nf = 100 (run BT4); displayed lines represent 
5, 15, 25, 35, 45, 55, 65, 75, 85, and 95 % of the shell mass. For the two runs, the upper frames 
show the corresponding effect of circularization of orbits; here the histograms record the number 
of fragments (crosses represent initial conditions, open circles the conditions at the end of the 
run) with orbit aspect ratio above Rmm/ Rmax, where i?min and i?max are consecutive minimum 
and maximum radii attained by a fragment along its orbit. 



of an average value of the Coulomb logarithm. We note that the fact that the value of 
the Coulomb logarithm thus determined remains of the same order of magnitude as that 
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Fig. 12. The energies associated with the system of minisatellites (left frame) and measurement 
of the Coulomb logarithm < ln{\) >— '^^.z'^^ ^i^W/^f (right frame) (run BS2)- 



obtained in Fig.nproves that the classical expression for the effects of dynamical friction 
captures the correct scaling with respect to the mass of the object subject to friction. 



6. Discussion and conclusions 

In this paper we have set up a numerical laboratory for the study of dynamical fric- 
tion in inhomogeneous quasi-spherical galaxies and made several experiments that have 
elucidated the physical mechanisms that participate in the process. 

The general character of the simulations presented here is distinctly different from 
that adopted in previous studies in this general research area, because we have addressed a 
physical scenario in which the slow collapse induced by dynamical friction approximately 
preserves the spherical symmetry present initially. This configuration is inspired by the 
astrophysical problem of studying the evolution of a spherical galaxy in the presence 
of a substantial globular cluster system or as a result of the capture of a large set of 
mini-satellites dragged in along random directions. 

From the theoretical point of view, the framework adopted here has an important 
advantage with respect to the traditional case of the study of the decay of the orbit of a 
single satellite, because it allows us to run smoother simulations that are basically free 
from other effects unrelated to dynamical friction, such as those associated with lack of 
equilibrium in the initial configuration. 

In the past, most of the attention has been focused on the effect of the host galaxy 
on the object subject to dynamical friction. Here, we have been able not only to test the 
adequacy of the classical formulae of dynamical friction (derived for ideal homogeneous 
models) in the inhomogeneous galaxy environment, but also to measure the evolution of 
the density and pressure tensor profiles in the galaxy induced by the presence of friction 
processes on a minority population of heavier particles. We have made quantitative tests 
that convince us that we are detecting genuine effects associated with secular evolution. 
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6.1. Plans for future investigations 

Starting from the encouraging quantitative results obtained in this paper, there are now 
two major issues that we would like to address in investigations planned for the near 
future, which appear to be approachable by the use of the numerical laboratory developed 
here: 

— In this paper wc have limited our study to the presence of one or two shells of heavy 
particles, so as to probe the reaction of the galaxy distribution function in a controlled 
way and to test whether the interaction involved can be described as "local" . The 
first immediate objective that we have in mind is to devise a set of more realistic sim- 
ulations, able to represent the distribution of the globular cluster system in galaxies 
better. 

— So far we have limited our discussion to the study of the orbital decay in an n = 3 
polytropic, isotropic galaxy model. Now we would like to test systematically the pro- 
cess of slow collapse induced by dynamical friction as a function of the concentration 
properties of the equilibrium model, for the case of more realistic galaxy models. 
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